m 



Fast Simulation Method 

for Parameter Reconstruction 

in Optical Metrology 

Sven Burger, ab Lin Zschiedrich, b Jan Pompom, b 
Frank Schmidt, ab Bernd Bodermann c 

Zuse Institute Berlin (ZIB), Takustrafie 7, D- 14 195 Berlin, Germany 

b JCMwave GmbH, Bolivarallee 22, D- 14 050 Berlin, Germany 

c Physikalisch-Teclmisclie Bundesanstalt Braunschweig (PTB), 
Bundesallee 100, D- 38 116 Braunschweig, Germany 

o' 

CN \ This paper will be published in Proc. SPIE Vol. 8681 (2013) 868119, (Metrology, Inspection, and Process 

?H ' Control for Microlithography XXVII, DOI: 10.1117/12.2011154), and is made available as an electronic preprint 

£~h| with permission of SPIE. One print or electronic copy may be made for personal use only. Systematic or multiple 

^h ■ reproduction, distribution to multiple locations via electronic or other means, duplication of any material in this 

\^0 | paper for a fee or for commercial purposes, or modification of the content of the paper are prohibited. 

^H ■ 

, ,! ABSTRACT 

A method for automatic computation of parameter derivatives of numerically computed light scattering signals 
is demonstrated. The finite-element based method is validated in a numerical convergence study, and it is 

©H" applied to investigate the sensitivity of a scatterometric setup with respect to geometrical parameters of the 
scattering target. The method can significantly improve numerical performance of design optimization, parameter 

tt . reconstruction, sensitivity analysis, and other applications. 
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1. INTRODUCTION 

t— I \ In optical metrology of nanostructures rigorous (i.e., accurate) simulation of light propagation is an essential 
^ > component PJ2l A challenge consists in reducing computation times for simulation results matching predefined 
accuracy requirements. This is especially important when real- world structures of complex geometry are consid- 
_~ ■ ered. 

We present a fast, finite-element based method to address such computation challenges. In this contribution 
t^J- \ we especially focus on finite-element based computation of derivatives of the propagating light fields (and of de- 
^^ ■ rived quantities like transmission or reflection intensities) with respect to geometrical parameters of the scattering 
target. As practical example we present a sensitivity analysis for patterns on a scatterometry reference stan- 
dard: dependence of the scatterometric signal on geometry parameters (CDs, sidewall-angles, corner-rounding) 
is evaluated in various parameter regimes. 

This paper is structured as follows: The background of our model is presented in Section [2j the numerical 
method is 
analysis o 
Section [5] 



5-H ■ method is described in Section [31 convergence results are reported in Section 2J and results of a sensitivity 
analysis of scatterometric signals from a pattern proposed as sample on a scatterometric standard is reported in 
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2. BACKGROUND / MODEL 

Light scattering off nanoscopic structures on scatterometry samples is modeled by the linear Maxwell's equations 
in frequency domain.™] From these a single equation for the electric field E can be derived: 



curl /i 1 curl E w 2 eE = iuiJ. 



(1) 



where e and \i arc the permittivity and permeability tensor, u is the time-harmonic frequency of the electromag- 
netic field, and the electric current J is source of an electromagnetic field. The domain of interest is separated 
into an infinite exterior J7 0Xt which hosts the given incident field and the scattered field, and an interior r2j nt where 
the total field is computed. Electromagnetic waves incident from the exterior to the interior at the boundaries 
between both domains are added to the right hand side of Eq. (|T|). For numerical simulations the infinite exterior 
is treated using transparent boundary conditions (using the perfectly matched layer method, PML). 

Transforming Eq. ([1} into weak formulation and discretizing it using finite elements yields a matrix equation: 



AE h 



f 



(2) 



where A is a sparse matrix, / contains the source terms, and E^ is the expansion of the electric field in a 
finite-dimensional FEM basis. 



Inversion of A and multiplication with the right hand side gives the solution E/ t : 

E„ = A- 1 / 



(3) 



Note that solutions corresponding to different sources incident on the same pattern can be obtained from the 
same inverted system matrix, given that A does not depend on the sources. E.g., when /i and ji correspond to 
incident light of two different polarizations, the corresponding near fields ~Eh,i and Eh,2 can be obtained from 
the same inverted system matrix: 



E/i : i 

E/j : 2 



= A-\ fl 
= A-\f 2 



(4) 
(5) 




Figure 1. Electric field intensity distribution in pseudo-color representation (red: high intensity, blue: low intensity). Top: 
linear color scale, bottom: logarithmic color scale. Left: S-polarized light, right: P-polarized light. 
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Figure 2. Schematics of the geometry of the investigated scatterometric target (unit cell of a ID-periodic grating). Free 
parameters of the model are the critical dimension (CD, width at h/2), the height h, pitch p, sidewall angle a, corner 
rounding radius R. 



Inversion of the system matrix (i.e., computation of A" 1 ) typically is the computationally most costly step, 
therefore re-using the same inverted matrix A^ 1 for TV sources reduces the computational costs approximately 
by a factor of A -1 , in a simulation setting where N independent source terms are present. 

In optimization problems, reconstruction problems and sensitivity studies, often an accurate measure of 
the partial derivative of the near field with respect to project parameters pi (e.g., geometry parameters, source 
parameters, material parameters), dpiEih, is required. As is well known, it is straight- forward in the finite-element 
context to compute these quantities by again re-using the inverted system matrix: 



5p,-E h = A- l [d Pi f - (d Pi A)E h ) 



Also higher-order derivatives 9™pjE/, can be computed, e.g., 



d 2 Pi E h = A- l [{d 2 pif) - (0 2 Pl A)E h - 2{dp i A){dp i E, h )) 



(6) 



(7) 



Here, d PiA is the Ath derivative of A with respect to parameter p,-, and d ptf is the Ath derivative of source 
term / with respect to parameter pi. 



3. NUMERICAL METHOD 

For rigorous simulations of the scattered light field we use the finite-element (FEM) Maxwell solver JCMsuite. 
This solver incorporates higher-order edge-elements, self-adaptive meshing, and fast solution algorithms for 
solving time-harmonic Maxwell's equations. Also, automatic computation of first- and higher-order parameter 
derivatives is implemented in the software. Previously the solver has, e.g., been used in scatterometric investiga- 
tions of EUV line masks (ID-periodic patterns), contact hole masks (2D-periodic patterns) and more complicated 
3D patterns P® Convergence studies in these investigations demonstrate that highly accurate, rigorous results 
can be attained even for the relatively large 3D computational domains which are typically present in 3D EUV 
setups. 

The workflow for the simulations is as follows: a scripting language (Matlab) automatically iterates the input 
parameter sets (physical parameters like geometrical dimensions and numerical parameters like mesh refinement). 
For each set, a triangular 2D mesh is created automatically by the built-in mesh generator. Then, the solver is 
started for computing the electromagnetic near field and its parameter derivatives, postprocessing is performed 
to extract, e.g., diffraction order efficiencies and their parameter derivatives, and results are evaluated and saved. 

Numerical settings which yield highly accurate results for the setup of interest in the presented investigations 
are identified in a convergence study (Section 2]). As numerical settings for the solver in the subsequent Section [3] 
on a sensitivity study, finite elements of third-order polynomial degree, and adaptive, error-estimator controlled 
meshing of the geometry in the computational domain and of transparent boundaries are chosen. This setting 
yields discrete problems with few ten thousands of unknowns (e.g., 30,000 unknowns), and few seconds (e.g., 
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Table 1. Parameter settings for the scatterometry standard simulations (compare Fig. \2§. Line height h, sidewall angle 
a, corner rounding radius R, illumination vacuum wavelength Ao, illumination inclination and rotation angle, 9 and <j>. 
Parameter settings for the scatterometry standard simulations (compare Fig. [2]) 




Figure 3. Finite-element mesh for spatial discretization of the geometry. Left: full geometry, right: detail at a rounded 
corner. 



4 sec) of computation time per computation (for computation of reflectivities and their parameter derivatives, 
for two polarizations, and for a specific physical setting). The FEM software solves these problems by direct LU 
factorization on a standard desktop computer. Figure Q] shows a graphical representation of a typical near-field 
intensity distribution. Please note that (as expected for this angle of incidence) the S-polarized incident wave 
leads to a smooth intensity distribution while the P-polarized incident wave leads to a highly discontinuous 
intensity distribution. 



4. MODEL VALIDATION 

In order to validate our model we perform a convergence study where we investigate how the computed quantities 
and their derivatives with respect to geometry parameters depend on the chosen numerical parameters. We 
investigate a geometry which could be used as part of a scatterometric standard.^ The investigated pattern is 
a ID-periodic line gratings etched into silicon (Si), with specific pitch (periodicity) and center line- width (CD) 
Figure [2] shows a schematics of the 2D setup for this test case. Table [T] shows parameter values of the project 
setup. Figure [3] shows a graphical representation of a 2D mesh. 

The pattern is illuminated from the superspace at oblique incidence with S- and P-polarized, monochromatic 
plane waves. The quantity of interest in this case is the intensity of light in the zero'th reflected diffraction order, 
Jo (Jo ~ |E| 2 , cf., Eq. QJ, and it's derivatives with respect to line width, height and sidewall angle, dlo/dcD, 
dlo/dh, dIo/d a , as function of varied geometry parameters. Please note that here, we normalize Jo with the 
intensity of the incoming light field, i.e., Iq is a dimensionless quantity. This numerical study is restricted to 
evaluation of intensities of the unpolarized light field, Jo, however, as the derivatives of the vectorial electric 
field amplitudes are computed (dpiEh, cf., Eq. [5]), also other quantities (sensitivities of all entries in the Miiller 
matrix) are accessible without extra computational costs. This numerical study is also restricted to ID-periodic 
patterns (i.e., 2D computational domains), however, the method can also be applied (and is implemented in the 
software) for 3D setups and/or isolated computational domains (i.e., non-periodic setups). 
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Figure 4. Dependence of the relative error of the reflectivity and its derivatives with respect to geometry parameters on 
numerical parameter p. Left: S-polarized incident light, Right: P-polarized incident light. 



Numerical errors as present in any numerical method for solving Maxwell's equations depend on the actual 
numerical settings. The two main numerical degrees of freedom for the finite-element method are the spatial 
discretization (mesh refinement) and the choice of ansatz functions which are used to approximate the fields on 
the spatial discretization mesh. The ansatz functions are typically defined by their polynomial degree p (when 
ansatz functions with a higher degree are chosen, this results in a larger basis for approximating the solution, and 
- more importantly - in higher approximation quality 3 ). Figure [5] shows how the numerical error of the reflection 
intensity and of its derivatives converges with finite element polynomial degree p. Relative errors are defined as 
normalized deviations from so-called quasi-exact results (results obtained at higher numerical discretization) ! 1 x '' 12 ' 
As can be seen from this Figure, very high levels of accuracy are reached for both, the reflected intensities, and 
for their derivatives with respect to geometry parameters. We have also checked that computing these derivatives 
using numerical differentiation yields the same numerical values (however, with worse convergence properties, 
and at significantly higher numerical cost). 



5. SENSITIVITY STUDY 

In order to demonstrate utility of the method we have performed several exemplary sensitivity studies. For 
given setups we investigate how the derivatives with respect to geometry parameters depend on specific physical 
parameter settings. This can be used to identify regimes where a scatterometric setup should work with higher 
sensitivity (yielding lower measurement uncertainties) than in other regimes. 

Figure [5] (left) shows how the scatterometric signal (zero order reflection intensity) varies with azimuthal 
angle of incidence of the illuminating plane waves. As expected, S- and P-polarization show a different behavior. 
The right part of this Figure shows how the sensitivity with respect to parameter variations depends on this 
angle. From this Figure it can, e.g., be seen that in this case sensitivity is about an order of magnitude higher 
for incident P-polarized light, and that absolute values of sensitivity are highest for small angles 9 (i.e., close to 
perpendicular incidence) . 

Figure [5] (left) shows how the scatterometric signal (zero order reflection intensity) varies with height of the 
grating lines. As in the previous case, S- and P-polarization show a different behavior. The right part of this 
Figure shows how the sensitivity with respect to parameter variations depends on the line height. From this 
Figure it can, e.g., be seen that in this case again, sensitivity is about an order of magnitude higher for incident 
P-polarized light, and that absolute values of sensitivity with respect to CD variations are highest for line of 
height h f=a 20 nm (in the investigated parameter regime) . 
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Figure 5. Left: Dependence of the scatterometric signal Jo on the azimuthal angle of incidence of the illuminating plane 
waves, for S- and P-polarization. Right: Dependence of the sensitivity with respect to parameter variations (CD, height, 
sidewall angle) on the angle of incidence. 
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Figure 6. Left: Dependence of the scatterometric signal Iq on the height of the grating lines h, for S- and P-polarization. 
Right: Dependence of the sensitivity with respect to parameter variations (CD, height, sidewall angle) on h. 



6. CONCLUSION 

To summarize, a method for automatic and computational-cost-effective computation of parameter derivatives of 
electromagnetic near fields and derived quantities has been demonstrated. This is useful for design optimization 
tasks, parameter reconstruction, sensitivity analysis, and other applications. A convergence study has been 
performed which demonstrated that very high levels of accuracy can be achieved. The method has been applied 
to investigate sensitivity of a scatterometric setup in different parameter regimes. 
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